Pythagoras goes Assembler

← Zurück zum Forum

Dr. Irata 10.01.2024 09:45

Wer kennt ihn nicht, den Satz des Pythagoras: "In allen rechtwinkligen Dreiecken ist die Summe der Flächen der Katheten- Quadrate gleich der Fläche des Quadrates der Hypotenuse"

A^2 + B^2 = C^2 ----> C = (A^2 + B^2)^0.5 also: C = Wurzel aus A Quadrat plus B Quadrat

Wir haben jetzt alles, um diese wichtige Formel in Assembler umzusetzen. Im Anschluß brauchen wie noch die p/q-Formel, sin/cos, e/ln

Nun zu Pythagoras (den brauche ich übrigens dringend für den Raycaster):

Das Problem ist in Assembler mit 8 Bit, daß man sehr schnell bei A^2 + B^2 auf eine 16-Bit Zahl kommt und dann auch noch die Wurzel dieser 16-Bit Zahl ziehen muss. Das Ergebnis kann wieder eine 8 Bit Zahl sein, aber zwischendurch muss man sehr wahrscheinlich durch die 16 Bit marschieren. Das ist jetzt aber alles gar kein Problem mehr. Wir bilden erst das Produkt AA, dann BB, macht dann eine 16-Bit Addiotion von beiden Produkten und davon dann ein 16-Bit Wurzel ziehen - und schwubs hat man das Ergebnis. Daraus bauen wir eine schöne kompakte Subroutine...

Liebe Grüße Peter

LarsImNetz 10.01.2024 13:19

Hi Peter, eigentlich suchst Du die Funktion hypot(a,b) Berechnen einer Hypotenuse. In meiner Compiler-Lib habe ich 2 Implementierungen gefunden, die Dir beide vielleicht nicht gefallen. Aber hier mal ein Beispiel:

; @see: https://arxiv.org/pdf/1904.09481.pdf

; hypot(a,b) returns sqrt(a^2+b^2) big problem is overrun
; if a == 0 then
;    h = b
;  else
;    r = b / a
;    h = a * sqrt(1 + r * r)
;  end if
; BUT: This works only with float values

Ich könnte meinen Code hier einfach abladen, aber: ich habe ein Problem mit den Parametern, wenn mein Compiler einen Funktionsaufruf macht:

 a := function_call(a,b,c)

Werden die Parameter a,b,c vom Compiler immer erst einmal auf 16 Bit "aufgeblasen" und dann auf einem eigenen Heap abgelegt: Dieser @heap_ptr liegt irgendwo im Speicher (Cassettenbuffer, oder am Ende des Sources...)

 lda a
 ldy #1
 sta (@heap_ptr),y
 lda #0
 iny
 sta (@heap_ptr),y
 lda b
 ldy #3
 sta (@heap_ptr),y
 lda b+1
 iny
 sta (@heap_ptr),y
 lda c
 ldy #5
 sta (@heap_ptr),y
 lda #0
 iny
 sta (@heap_ptr),y
 jsr function_call_iii
 ...

In allen selbst geschriebenen Funktionen fängt mein Source deshalb meist erstmal an die Parameter auszulesen:

function_call_iii
 ldy #1
 lda (@heap_ptr),y
 sta parameter_a
 ldy #3
 lda (@heap_ptr),y
 sta parameter_b
 iny
 lda (@heap_ptr),y
 sta parameter_b+1
 iny
 lda (@heap_ptr),y
 sta parameter_c 
 ...

Das ist straightforward für meinen Compiler, diesen ganzen .proc krams (von Mads/cc65) kann der atasm Assembler nicht, den ich ja vorzugsweise verwende. Des Weiteren kommen Rückgabewerte immer in den X-Y-Registern zurück, immer der Form: (y+x*256)

Desweiteren gibt es 'lokale' Variablen, diese beginnen immer mit '?'. Und ein '.local' ändert intern im Atasm nur einen Zähler, so das Variablen, die mit '?' beginnen solange gleich sind, bis ein erneutes '.local' kommt.

Soviel zu meinen Routinen.

Hier eine Implementierung:

 .local
 ?POW2 = @REG+1

; Funktion berechnet die Hypotenuse von 2 Parametern (jeweils 8 Bit gross)
@HYPOT
@HYPOT_II
  LDY #1
  LDA (@HEAP_PTR),Y ; wert a holen
  TAX

  LDY #3
  LDA (@HEAP_PTR),Y ; wert b holen
  TAY

  clc
  lda @pow2_low,x
  adc @pow2_low,y
  sta ?pow2
  lda @pow2_high,x
  adc @pow2_high,y
  sta ?pow2+1
; TODO: wenn carry gesetzt ist, ist der Wert > 65535
;  lda #0
;  adc #0
;  sta ?pow2+2

  jmp @SQRT_INNER

@pow2_low
 .byte <0, <1, <4, <9, <16, <25, <36, <49, <64, <81, <100, <121, <144, <169, <196, <225
 .byte <256, <289, <324, <361, <400, <441, <484, <529, <576, <625, <676, <729, <784, <841, <900, <961
 .byte <1024, <1089, <1156, <1225, <1296, <1369, <1444, <1521, <1600, <1681, <1764, <1849, <1936, <2025, <2116, <2209
 .byte <2304, <2401, <2500, <2601, <2704, <2809, <2916, <3025, <3136, <3249, <3364, <3481, <3600, <3721, <3844, <3969
 .byte <4096, <4225, <4356, <4489, <4624, <4761, <4900, <5041, <5184, <5329, <5476, <5625, <5776, <5929, <6084, <6241
 .byte <6400, <6561, <6724, <6889, <7056, <7225, <7396, <7569, <7744, <7921, <8100, <8281, <8464, <8649, <8836, <9025
 .byte <9216, <9409, <9604, <9801, <10000, <10201, <10404, <10609, <10816, <11025, <11236, <11449, <11664, <11881, <12100, <12321
 .byte <12544, <12769, <12996, <13225, <13456, <13689, <13924, <14161, <14400, <14641, <14884, <15129, <15376, <15625, <15876, <16129
 .byte <16384, <16641, <16900, <17161, <17424, <17689, <17956, <18225, <18496, <18769, <19044, <19321, <19600, <19881, <20164, <20449
 .byte <20736, <21025, <21316, <21609, <21904, <22201, <22500, <22801, <23104, <23409, <23716, <24025, <24336, <24649, <24964, <25281
 .byte <25600, <25921, <26244, <26569, <26896, <27225, <27556, <27889, <28224, <28561, <28900, <29241, <29584, <29929, <30276, <30625
 .byte <30976, <31329, <31684, <32041, <32400, <32761, <33124, <33489, <33856, <34225, <34596, <34969, <35344, <35721, <36100, <36481
 .byte <36864, <37249, <37636, <38025, <38416, <38809, <39204, <39601, <40000, <40401, <40804, <41209, <41616, <42025, <42436, <42849
 .byte <43264, <43681, <44100, <44521, <44944, <45369, <45796, <46225, <46656, <47089, <47524, <47961, <48400, <48841, <49284, <49729
 .byte <50176, <50625, <51076, <51529, <51984, <52441, <52900, <53361, <53824, <54289, <54756, <55225, <55696, <56169, <56644, <57121
 .byte <57600, <58081, <58564, <59049, <59536, <60025, <60516, <61009, <61504, <62001, <62500, <63001, <63504, <64009, <64516, <65025

@pow2_high
 .byte >0, >1, >4, >9, >16, >25, >36, >49, >64, >81, >100, >121, >144, >169, >196, >225
 .byte >256, >289, >324, >361, >400, >441, >484, >529, >576, >625, >676, >729, >784, >841, >900, >961
 .byte >1024, >1089, >1156, >1225, >1296, >1369, >1444, >1521, >1600, >1681, >1764, >1849, >1936, >2025, >2116, >2209
 .byte >2304, >2401, >2500, >2601, >2704, >2809, >2916, >3025, >3136, >3249, >3364, >3481, >3600, >3721, >3844, >3969
 .byte >4096, >4225, >4356, >4489, >4624, >4761, >4900, >5041, >5184, >5329, >5476, >5625, >5776, >5929, >6084, >6241
 .byte >6400, >6561, >6724, >6889, >7056, >7225, >7396, >7569, >7744, >7921, >8100, >8281, >8464, >8649, >8836, >9025
 .byte >9216, >9409, >9604, >9801, >10000, >10201, >10404, >10609, >10816, >11025, >11236, >11449, >11664, >11881, >12100, >12321
 .byte >12544, >12769, >12996, >13225, >13456, >13689, >13924, >14161, >14400, >14641, >14884, >15129, >15376, >15625, >15876, >16129
 .byte >16384, >16641, >16900, >17161, >17424, >17689, >17956, >18225, >18496, >18769, >19044, >19321, >19600, >19881, >20164, >20449
 .byte >20736, >21025, >21316, >21609, >21904, >22201, >22500, >22801, >23104, >23409, >23716, >24025, >24336, >24649, >24964, >25281
 .byte >25600, >25921, >26244, >26569, >26896, >27225, >27556, >27889, >28224, >28561, >28900, >29241, >29584, >29929, >30276, >30625
 .byte >30976, >31329, >31684, >32041, >32400, >32761, >33124, >33489, >33856, >34225, >34596, >34969, >35344, >35721, >36100, >36481
 .byte >36864, >37249, >37636, >38025, >38416, >38809, >39204, >39601, >40000, >40401, >40804, >41209, >41616, >42025, >42436, >42849
 .byte >43264, >43681, >44100, >44521, >44944, >45369, >45796, >46225, >46656, >47089, >47524, >47961, >48400, >48841, >49284, >49729
 .byte >50176, >50625, >51076, >51529, >51984, >52441, >52900, >53361, >53824, >54289, >54756, >55225, >55696, >56169, >56644, >57121
 .byte >57600, >58081, >58564, >59049, >59536, >60025, >60516, >61009, >61504, >62001, >62500, >63001, >63504, >64009, >64516, >65025

Aufmerksame Leser sehen sofort, das SQRT() fehlt, das kommt hier:

 .local
; fast integer SQRT

?THI = @REG
?MLO = @REG+1
?MHI = @REG+2

; find more here https://www.codebase64.org/doku.php?id=base:fast_sqrt
@SQRT
@SQRT_I
      LDY #1
      LDA (@HEAP_PTR),Y
      STA ?MLO
      INY
      LDA (@HEAP_PTR),Y
      STA ?MHI

        CLC         ; clear bit 16 of M
@SQRT_INNER
        LDY #$00    ; R = 0
        LDX #$07
?loop
        TYA
        ORA ?stab-1,X
        STA ?THI     ; (R ASL 8) | (D ASL 7)
        LDA ?MHI
        BCS ?skip0  ; M >= 65536? then T <= M is always true
        CMP ?THI
        BCC ?skip1  ; T <= M
?skip0
        SBC ?THI
        STA ?MHI     ; M = M - T
        TYA
        ORA ?stab,x
        TAY         ; R = R OR D
?skip1
        ASL ?MLO
        ROL ?MHI     ; M = M ASL 1
        DEX
        BNE ?loop

        ; last iteration

        BCS ?skip2
        STY ?THI
        LDA ?MLO
        CMP #$80
        LDA ?MHI
        SBC ?THI
        BCC ?skip3
?skip2
        INY         ; R = R OR D (D is 1 here)
?skip3
        RTS

?stab   .BYTE $01,$02,$04,$08,$10,$20,$40,$80

LG Lars

LarsImNetz 10.01.2024 13:55

BTW: Ein Test meinerseits hat ergeben, das die hier angegebene Routine für die Wurzel ca. 45% schneller ist als die von Lee Davison. Aber probiere es selbst aus. Ich habe es auch nicht geglaubt.

Dr. Irata 10.01.2024 13:57

Hallo Lars, ja so etwas suchte ich - interessanter Ansatz muss ich sagen, sieht aber recht kompliziert aus. Allerdings haben wir jetzt alles, was man braucht, um c stabil auszurechnen - eine schnelle Multiplikationsroutine mit 16-Bit Ergebnis, die Addition zweier 16-Bit-Zahlen habe ich vorhin gemacht (das ist recht einfach und kurz) und dann diese schnelle Wurzelfunktion einer 16 Bit Zahl. Damit kann man nun recht einfach c ausrechnen!! Ich stelle das in Kürze hier rein, dann können wir zusammen den Code noch optimieren und haben es für gut! Liebe Grüße Peter

p.s. An gute und schnelle Algorithmen bin ich immer sehr interessiert!!

Dr. Irata 10.01.2024 13:58

... 45% schneller ist natürlich mega interessant!!

Die Berechnung der Hypothenuse mit r=a/b und dann Wurzel aus (1+r^2) ist ja nur umgeformter Pythagoras. Der Vorteil dabei ist, daß man mit deutlich niedrigeren Zahlen rechnen kann, da man nun nur noch das Quadrat eines Quotienten bilden muss. a=20 und b=10 würde normal Wurzel aus (400+100) also Wurzel aus (500) = 22 Rest 16 = 22,36 Bildet man vorher r=a/b = 20/10 = 2 und dann Wurzel (1+ r*r) also Wurzel 5 = 2,236 * 10 = 22,36 hat man deutlich kleinere Zahlen, allerdings wird man wenn man kein Fließkomma nutzt einen höheren Fehler produzieren!

Dr. Irata 10.01.2024 14:00

... kannst du die Routine noch in reinen Assembler übersetzten? Ich verstehe die Syntax mit dem ? und dem @ nicht...

LarsImNetz 10.01.2024 15:20

Das ist reiner MAC65 Code, das @ hat keinerlei extra Bedeutung. Es ist wie ein Buchstabe. Auch MADS erlaubt es als ersten Wert in einer Variablen, aber hier nochmal der komplette Code, etwas mehr kommentiert.

.local

REG=214        ; irgendwo in der Zeropage
_HEAP_PTR=1024 ; wo etwas Platz ist im Speicher

POW2=REG+1    ; + 1 weil SQRT dort den Wert erwartet zum Ziehen der Wurzel

; Funktion berechnet die Hypotenuse von 2 Parametern (jeweils 8 Bit gross)
HYPOT
HYPOT_II
  LDY #1
  LDA (_HEAP_PTR),Y ; wert a holen
  TAX

  LDY #3
  LDA (_HEAP_PTR),Y ; wert b holen
  TAY

; Berechnung von n*n mittels Tabelle, schneller geht es nicht, kostet aber 512 Bytes Speicher!
  CLC
  LDA pow2_low,x
  ADC pow2_low,y
  STA pow2
  LDA pow2_high,x
  ADC pow2_high,y
  STA pow2+1     ; pow2 ist x*x + y*y


; TODO: wenn carry gesetzt ist, ist der Wert > 65535
;  lda #0
;  adc #0
;  sta pow2+2

  JMP SQRT_INNER

; Enthaelt für den Index n das Low-Byte von n*n 
pow2_low
 .byte <0, <1, <4, <9, <16, <25, <36, <49, <64, <81, <100, <121, <144, <169, <196, <225
 .byte <256, <289, <324, <361, <400, <441, <484, <529, <576, <625, <676, <729, <784, <841, <900, <961
 .byte <1024, <1089, <1156, <1225, <1296, <1369, <1444, <1521, <1600, <1681, <1764, <1849, <1936, <2025, <2116, <2209
 .byte <2304, <2401, <2500, <2601, <2704, <2809, <2916, <3025, <3136, <3249, <3364, <3481, <3600, <3721, <3844, <3969
 .byte <4096, <4225, <4356, <4489, <4624, <4761, <4900, <5041, <5184, <5329, <5476, <5625, <5776, <5929, <6084, <6241
 .byte <6400, <6561, <6724, <6889, <7056, <7225, <7396, <7569, <7744, <7921, <8100, <8281, <8464, <8649, <8836, <9025
 .byte <9216, <9409, <9604, <9801, <10000, <10201, <10404, <10609, <10816, <11025, <11236, <11449, <11664, <11881, <12100, <12321
 .byte <12544, <12769, <12996, <13225, <13456, <13689, <13924, <14161, <14400, <14641, <14884, <15129, <15376, <15625, <15876, <16129
 .byte <16384, <16641, <16900, <17161, <17424, <17689, <17956, <18225, <18496, <18769, <19044, <19321, <19600, <19881, <20164, <20449
 .byte <20736, <21025, <21316, <21609, <21904, <22201, <22500, <22801, <23104, <23409, <23716, <24025, <24336, <24649, <24964, <25281
 .byte <25600, <25921, <26244, <26569, <26896, <27225, <27556, <27889, <28224, <28561, <28900, <29241, <29584, <29929, <30276, <30625
 .byte <30976, <31329, <31684, <32041, <32400, <32761, <33124, <33489, <33856, <34225, <34596, <34969, <35344, <35721, <36100, <36481
 .byte <36864, <37249, <37636, <38025, <38416, <38809, <39204, <39601, <40000, <40401, <40804, <41209, <41616, <42025, <42436, <42849
 .byte <43264, <43681, <44100, <44521, <44944, <45369, <45796, <46225, <46656, <47089, <47524, <47961, <48400, <48841, <49284, <49729
 .byte <50176, <50625, <51076, <51529, <51984, <52441, <52900, <53361, <53824, <54289, <54756, <55225, <55696, <56169, <56644, <57121
 .byte <57600, <58081, <58564, <59049, <59536, <60025, <60516, <61009, <61504, <62001, <62500, <63001, <63504, <64009, <64516, <65025

; Enthaelt für den Index n das High-Byte von n*n 
pow2_high
 .byte >0, >1, >4, >9, >16, >25, >36, >49, >64, >81, >100, >121, >144, >169, >196, >225
 .byte >256, >289, >324, >361, >400, >441, >484, >529, >576, >625, >676, >729, >784, >841, >900, >961
 .byte >1024, >1089, >1156, >1225, >1296, >1369, >1444, >1521, >1600, >1681, >1764, >1849, >1936, >2025, >2116, >2209
 .byte >2304, >2401, >2500, >2601, >2704, >2809, >2916, >3025, >3136, >3249, >3364, >3481, >3600, >3721, >3844, >3969
 .byte >4096, >4225, >4356, >4489, >4624, >4761, >4900, >5041, >5184, >5329, >5476, >5625, >5776, >5929, >6084, >6241
 .byte >6400, >6561, >6724, >6889, >7056, >7225, >7396, >7569, >7744, >7921, >8100, >8281, >8464, >8649, >8836, >9025
 .byte >9216, >9409, >9604, >9801, >10000, >10201, >10404, >10609, >10816, >11025, >11236, >11449, >11664, >11881, >12100, >12321
 .byte >12544, >12769, >12996, >13225, >13456, >13689, >13924, >14161, >14400, >14641, >14884, >15129, >15376, >15625, >15876, >16129
 .byte >16384, >16641, >16900, >17161, >17424, >17689, >17956, >18225, >18496, >18769, >19044, >19321, >19600, >19881, >20164, >20449
 .byte >20736, >21025, >21316, >21609, >21904, >22201, >22500, >22801, >23104, >23409, >23716, >24025, >24336, >24649, >24964, >25281
 .byte >25600, >25921, >26244, >26569, >26896, >27225, >27556, >27889, >28224, >28561, >28900, >29241, >29584, >29929, >30276, >30625
 .byte >30976, >31329, >31684, >32041, >32400, >32761, >33124, >33489, >33856, >34225, >34596, >34969, >35344, >35721, >36100, >36481
 .byte >36864, >37249, >37636, >38025, >38416, >38809, >39204, >39601, >40000, >40401, >40804, >41209, >41616, >42025, >42436, >42849
 .byte >43264, >43681, >44100, >44521, >44944, >45369, >45796, >46225, >46656, >47089, >47524, >47961, >48400, >48841, >49284, >49729
 .byte >50176, >50625, >51076, >51529, >51984, >52441, >52900, >53361, >53824, >54289, >54756, >55225, >55696, >56169, >56644, >57121
 .byte >57600, >58081, >58564, >59049, >59536, >60025, >60516, >61009, >61504, >62001, >62500, >63001, >63504, >64009, >64516, >65025

 .local

; fast integer SQRT

THI = REG
MLO = REG+1

; find more here https://www.codebase64.org/doku.php?id=base:fast_sqrt
SQRT
SQRT_I
      LDY #1
      LDA (_HEAP_PTR),Y
      STA MLO
      INY
      LDA (_HEAP_PTR),Y
      STA MLO+1

      CLC         ; clear bit 16 of M
SQRT_INNER
      LDY #$00    ; R = 0
      LDX #$07
loop
      TYA
      ORA stab-1,X
      STA THI     ; (R ASL 8) | (D ASL 7)
      LDA MLO+1
      BCS skip0  ; M >= 65536? then T <= M is always true
      CMP THI
      BCC skip1  ; T <= M
skip0
      SBC THI
      STA MLO+1     ; M = M - T
      TYA
      ORA stab,x
      TAY         ; R = R OR D
skip1
      ASL MLO
      ROL MLO+1     ; M = M ASL 1
      DEX
      BNE loop

      ; last iteration

      BCS skip2
      STY THI
      LDA MLO
      CMP #$80
      LDA MLO+1
      SBC THI
      BCC skip3
skip2
      INY         ; R = R OR D (D is 1 here)
skip3
; Das Ergebnis ist im Register (y+x*256)
      RTS

stab
      .BYTE $01,$02,$04,$08,$10,$20,$40,$80
Dr. Irata 11.01.2024 15:26

... und hier nun meine Routine für die Berechnung von c nach Pythagoras --> Wurzel aus (aa+bb) :

        ; Pythagoras c = (a*a + b*b)^0,5
        ; a=var1
        ; b=var2
        ; c=var3

        .proc pythagoras
        ; a^2 ausrechnen
        ;............... var1 * var1 = temp1/temp1+1 .......
        ;.... Multiplikation 8 Bit * 8 Bit auf 16 Bit ....
        lda var1
        sta temp1
        lda #0
        ldx #8
        lsr temp1
lp1
        bcc lp2
        clc
        adc var1
lp2
        ror
        ror temp1
        dex
        bne lp1
        sta temp1+1

        ; b^2 ausrechnen
        ;............... var2 * var2 = temp2/temp2+1 .......
        ;.... Multiplikation 8 Bit * 8 Bit auf 16 Bit ....
        lda var2
        sta temp2
        lda #0
        ldx #8
        lsr temp2
lp3
        bcc lp4
        clc
        adc var2
lp4
        ror
        ror temp2
        dex
        bne lp3
        sta temp2+1

        ; 16 Bit Addition
        ; temp1/temp1+1  +  temp2/temp2+1  =  temp3/temp3+1
        clc
        lda temp1
        adc temp2
        sta temp3
        lda temp1+1
        adc temp2+1
        adc #0
        sta temp3+1

        ; Quadratwurzel 16 Bit
        ; nach Lee Davison
        ; Quadratwurzel temp3/temp3+1 ---> Ergebnis in var3
        ldx #8
lp5
        sec
        lda temp3+1
        sbc #64
        tay
        lda temp4
        sbc var3
        bcc lp6
        sty temp3+1
        sta temp4
lp6
        rol var3
        asl temp3
        rol temp3+1
        rol temp4
        asl temp3
        rol temp3+1
        rol temp4
        dex
        bne lp5
        rts
        .endp
        ;das Ergebnis steht in var3

... viel Spass beim Testen und Optimieren!!

Liebe Grüße Peter

Dr. Irata 02.03.2024 08:47

Die Routine ist ja als Subroutine "gekapselt" - das kann man zum testen natürlich auch anders machen und das .proc / .endp weglassen. Außerdem müssen natürlich alle Variablen vorher initialisiert werden im Assembler z.B.
var1 .byte 00